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Abstract 

The properties of a dilute granular gas in the homogeneous cooling state are mapped to those of a 
stationary state by means of a change in the time scale that does not involve any internal property 
of the system. The new representation is closely related with a general property of the granular 
temperature in the long time limit. The physical and practical implications of the mapping are 
discussed. In particular, simulation results obtained by the direct simulation Monte Carlo method 
applied to the scaled dynamics are reported. This includes ensemble averages and also the velocity 
autocorrelation function, as well as the self-diffusion coefficient obtained from the latter by means 
of the Green-Kubo representation. In all cases, the obtained results are compared with theoretical 
predictions. 

PACS numbers: PACS Numbers: 45.70.-n,51.10.+y,05.20.Dd 
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I. INTRODUCTION 



A granular fluid is a collection of macroscopic particles interacting via short range hard 
inelastic collisions. Particles move in a ballistic way between collisions and total momentum 

n 

is conserved [1] . The prototypical idealized model for granular fluids is a system of inelastic 
smooth hard spheres or disks, with the inelasticity of collisions being described by means 
of a constant (independent of the relative velocity) coefficient of normal restitution. Then, 
in the last years the traditional methods of kinetic theory and non-equilibrium statistical 
mechanics have been extended to the case of inelastic collisions. Quite remarkably, it has 
been realized that the single feature of incorporating energy dissipation in collisions is able to 
provide a theoretical scheme where many of the peculiar features exhibited by real granular 
fluids can be tackled. This includes phenomena such as the development of strong density 
and temperature inhomogeneities that are not induced by the boundary conditions 3, Q|, 
spontaneous symmetry breaking in partitioned jj, |^ and non-partitioned systems JgI], seg- 
regation in systems composed of different kind of particles [3], and pattern formation 8], to 
cite a few examples. In most of these cases, the usefulness of a collective description of the 
system in terms of hydrodynamic-like equations has been verified. Such a description can 
only be fully understood and justified by starting from a more fundamental particle level, 
as considered in kinetic theory. 

Due to energy dissipation in collisions, granular systems do not present a stationary, 
homogeneous and isotropic state similar to the equilibrium one of ordinary fiuids. The 
simplest possible state corresponds to a freely evolving homogenous and isotropic system 
whose energy decays monotonically in time, the so-called homogeneous cooling state (HCS). 
This state plays a relevant role in order to investigate the transport properties of a granular 
fiuid, since it provides the zeroth order in the gradients approximation when applying the 
Chapman-Enskog procedure to derive hydrodynamics from a given kinetic equation 



ing the 

HQ- 



Also, linear response around this state has been studied and formal expressions for the 
Navie.Sto.es t^a^spo. coeffle.en. Mve been denve. HQ. T.ey are the generalization 
to inelastic systems of the well-known Green-Kubo formulas. Besides their theoretical inter- 
est, they allow a direct determination of the transport coefficients from the dynamics of the 
system in the HCS, without introducing any additional approximation, by using numerical 
simulation methods. 
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Molecular dynamics (MD) simulation provides a method to investigate a system of parti- 
cles at the most fundamental level of description. Nevertheless, when applied to a granular 
fluid in the HCS, several limitations show up. First, since the system is continuously cooling, 
the typical velocity of the particles becomes very small rather soon and numerical inaccu- 
racies become very large. In principle, this could be solved by introducing some kind of 
external thermostat, but then the relationship between the original HCS and the state be- 
ing actually simulated is not clear. Another possibility is to take advantage of the fact that 
there is no intrinsic time scale in a system of hard particles, and to rescale the velocity of all 
particles after every collision, so that the energy is forced to remain constant (l^. Although 
it seems that this method must lead to correct results for time-independent properties of 
the HCS, e.g. structural properties or the own scaled velocity distribution, in the infinite 
system limit, it is not evident how to extract from the simulation data properties of the 
actual dynamics of the system involving time fluctuations or two-time correlations. 

Very recently, a procedure has been introduced according to which the dynamics of the 

system in the HCS is exactly mapped onto the dynamics around a steady state by means 

I I 

of a change in the time scale being used 1.14] • The change is independent of the state of 
the system. This is possible because the temperature of the HCS becomes independent of 
its initial value in the long time limit. This is a very strong and fundamental property of 
that state that has not received too much attention up to now. Then, the existence of the 
steady state in the scaled time representation is tied to the own physical properties of the 
mechanism of energy dissipation. 

The second limitation of the MD simulation of a granular fluid in the HCS is associated 
with the fact that this state is unstable with respect to spacial long wavelength perturba- 
tions . This instability has been identifled in the context of hydrodynamics, and an 
expression for the critical size of the system, beyond which it becomes unstable, has been 
obtained. The critical size is a function of the density and the coefficient of restitution, 
decreasing with the former and increasing with the latter. In practice, this implies that for 
high densities and/or small values of the coefficient of restitution, only very small systems 
can be simulated in the HCS, and undesired finite size effects might influence the results. 

The aim of this paper is to investigate in detail the physical and practical implications 
of the steady state representation of the HCS mentioned above for a low density granular 
gas. Attention will be paid not only to the one-time properties of the system, but also 
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to two-time correlation functions. While the former can be discussed on the basis of the 
inelastic Boltzmann equation, the analysis of the correlation functions requires to introduce 
an equation for the two-particle and two-time distribution function. This is done b y a direct 
extension of the methods used in the elastic case for out of equilibrium systems 16|. It is 
shown that both kind of properties can be expressed in terms of averages over the stationary 
state of the system. Special emphasis will be put on the relationship between the theoretical 
description of the system in terms of reduced distribution functions in the low density 
limit and the underlying A^-particle dynamics. This is important in order to implement 
the calculation of a given property by means of the direct simulation Monte Carlo method 
(DSMC) j^. It must be kept in mind that this method is designed not just as a numerical 
tool to solve the Boltzmann equation, but as a real A^-particle dynamics simulation of a low 
density gas. In this sense, it is expected to provide not only the one-particle distribution 
function of the system, but its complete dynamical description. 

One of the main practical advantages of the DSMC method is that it allows to incorpo- 
rate at the level of the particle dynamics the spatial symmetry properties of the particular 
physical situation of interest. For instance, it is easy to restrict the dynamics of the system 
so that it remains homogeneous, with no possibility of developing spacial inhomogeneities. 
In this way, the second limitation of MD simulations of the HCS can be overcome for a low 
density gas with no restriction on the size of the system or the number of particles being 
used. The combination of this feature of the DSMC method and the steady representation 
of the HCS offers an almost unique way to investigate the properties of a dilute granular 
gas. 

The plan of the paper is as follows. In Sec. m the Boltzmann equation and some of 
the results for the one-particle distribution function of the HCS are shortly reviewed. The 
peculiar long time behavior of the temperature is indicated, and in Sec. IIIII it is used to 
construct a representation of the dynamics of the system in which the distribution function 
of the HCS becomes time independent. It is shown that the asymptotic steady temperature is 
determined by an intrinsic property of the system, namely the cooling rate. Then, dynamical 
simulations of a system of inelastic hard disks with the DSMC method are presented, and 
the numerical results for the cooling rate obtained from the values of the steady temperature 
are compared with the existing theoretical predictions. Results for the fourth moment of 
the velocity distribution function are also reported. They are consistent with those obtained 
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previously by using the original dynamics, in which the HCS is time dependent. 

Time correlation functions of dynamical variables in the HCS are addressed in Sec. IIV| 
where their low density limit is analyzed and the relationship between their values in the 
original and scaled dynamics is established. The particular case of the correlation functions 
appearing in the Green-Kubo form of the transport coefficients of a low density gas 
is considered. As an example, the velocity autocorrelation function, and from it the self- 
diffusion coefficient, are computed in Sec. |3for a system of inelastic hard disks, and the 
results are compared with the theoretical predictions obtained by the Chapman- Enskog 
method in the first Sonine approximation. The paper ends with a summary and a brief 
discussion. 



II. BASIC EQUATIONS AND THE HOMOGENEOUS COOLING STATE 

We consider a system of inelastic hard spheres {d = 3) or disks {d = 2) of mass m and 
diameter a. The position and velocity coordinates of particle i at time t will be represented 
by Ri{t) and Vi{t), respectively. The particle dynamics consists of free streaming until a 
pair of particles i and j are at contact and a collision takes place. The effect of the collision 
is to instantaneously change the velocities of the two involved particles according to 

v,,^v: = b^v, = V, - (a ■ V,,) a 

V, - v; ^ h^v, = v, + ^{a. a, (i) 

where Vij = Vi — Vj is the relative velocity, a is the unit vector pointing from the center of 
particle j to the center of particle i at contact, and a is the coefficient of normal restitution. 
It is defined in the range < a < 1 and will be considered as velocity independent along 
this work. 

In the low density limit, it is assumed that the time evolution of the one-particle distri- 
bution function of the system f{r,v,t) is accurately described by the inelastic Boltzmann 



equation 
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\dt dr ^ 

with the inelastic Boltzmann collision operator given by 



^ + v-4-]fir,v,t) = J[f,f], (2) 



J[f, f] ^ a'-' J dv, J da e{g ■a)g a (a-%-' - l) /(r, v, t)f{r, v,, t). (3) 



Here g = v — Vi, is the Heaviside step function, and b^^{v,Vi) is an operator replacing 
all the velocities v and vi appearing to its right by the precoUisional values v* and vl given 
by 

,-1 l + a 
V = b„v = v cr ■ Qlcr, 

vl = b^'vi = v, + ^^{d- ■ g)a. (4) 

We are using lower-case symbols to represent field variables as those appearing in the one- 
particle distribution function, while capital symbols are used for the particle variables. Equa- 
tion has a particular solution describing the HCS and having the scaling property 

fncsiv, t) = nHVQ^{t)xHcs{c), (5) 

where 



vo{t) 



2kBTHcs{ty'^' 



m 

is the thermal velocity and Xhcs{<^) is an isotropic function of 

V 



(6) 



(7) 



vo{t) ■ 

In Eq. (jni), ks is the Boltzmann constant that is usually set equal to unity in the literature 
of granular fiows. We prefer to keep it here for dimensional reasons. The number of particles 
density n and the granular temperature T(t) are defined, for arbitrary situations, in terms 
of the one-particle distribution function in the usual way. 



n{r,t) = J dvf{r,v,t), (8) 

^n{r, t)kBT{r, t) = j dv^-m{v - uff{r, v, t), (9) 

n{r,t)u{r,t) = J dvvf{r,v,t). (10) 

Therefore, in the HCS all the time dependence of the distribution function occurs through 
the temperature Tucsit), whose evolution equation is easily obtained from the own Boltz- 
mann equation 

^Tncsit) + CHCs{t)THcs{t) = 0, (11) 
where the cooling rate (hcs is given by 

(Hcsit) = ^ {d±3\/ I '^''1 "^""^ " c^\'xhcs{c)xhcs{ci). (12) 

2r ( — ) a 



Then, (hcs is proportional to T^/^g and, therefore, Eq. (fTT|) can be formally integrated to 
get the explicit time dependence of the temperature in the HCS, 

-2 



Tncsit) = Thcs{0) [l + ^CHCsiO)t^ ■ (13) 



This expression is known as Haff 's law |20|] and it has the interesting property of becoming 
independent of the initial temperature in the long time limit. More precisely, it reduces to 

Tncs-^mr', (14) 

with 



/,N /)ml/2 /.N ^ ' 



In the last transformation we have introduced the time-independent dimensionless cooling 
rate 

Co = (16) 

where i = {nncr'^'^)^^ is proportional to the mean free path. The above long time behavior 
of the temperature implies the existence of an asymptotic regime in which the HCS becomes 
independent of the initial condition or, in other words, all the homogeneous cooling states of 
a given system tend to converge in the long time limit. This is an exact property following 
from the own existence of the HCS and it will be exploited in the next Section in order to 
introduce a steady representation of the HCS. 

Substitution of Eq. (0) into the Boltzmann equation Q provides a closed integro- 
differential equation for the function Xhcs{c), 

Co 9 

'idc ' ^^^^^^^ " Jc[xhcs, Xhcs], (17) 

Jc[Xhcs,Xhcs\ = j ^'^^ J daQ[{c-ci)-a]{c-ci)-a[a'%:^^{c,Ci) - l\ Xhcs{c)xhcs{ci). 

(18) 

The operator h^^{c, Ci) is again defined by Eqs. (jH) but substituting the velocities Vi by 
c, Ci. 

The solution of Eq. p7|) is only partially known. In particular, an expansion in Sonine 
polynomials has been considered. To first order, Xhcs{.c) is approximated by 



Xhcs[c) - L 



l + a2(«)^(')(c')], (19) 



where 



2^ _ d + 2 2 d{d + 2) 
" ~2 ^ 8 ' 



(20) 



The coefficient 02(0;) turns out to be proportional to the fourth cumulant of the distribution, 
namely 

= ^[(eV^^l. (21) 



d{d + 2 

{c') ^ J dcc'xHcsic). (22) 
When Eq. (fTI?|) is substituted into Eq. (fTTj) a closed equation for 02 is obtained. If nonlinear 



terms in 02 are neglected in this equation, it is found 



2l|: 



a2{a] 



16(l-a)(l-2a2 



9 + 2Ad+ {8d - 41) + SOa^ _ ^q^s 
In the same approximation, Eqs. ()12|1 and pfi|) yield 



Co 



1 + -77:^2(0) 

lb 



(23) 



(24) 



The above expression for Xhcs{c) is expected to be accurate in the thermal velocity 
region, i.e. for velocities c of the order of a few units. This has been conffimed by DSMC 



simulations of a granular gas 



221. 



III. STEADY STATE REPRESENTATION OF THE HCS 



The long time behavior of the HCS discussed in the preceding Section, suggests that a 
steady representation of it can be obtained on a time scale r defined by 



ujqt = In — 
to 



(25) 



where cuq and to are arbitrary constants. Consistently, the velocity Wi of particle i is given 
in the new scale by 

W,(r) = cuotoe"''"^(t) = uJom{t). (26) 

The particle dynamics in these variables consists of an accelerating streaming between col- 
lisions, 



d 

d' 
d_ 



— i?,(r) = W.ir), 



(27) 



while the effect of a colhsion between particles i and j is to instantaneously modify their 
velocities accordingly with the same rules as given in Eqs. (Q), i.e. 

1 + a 
\ + a 

Wj ^ w; = b^Wj = Wj + ^— (a ■ w^j) a. (28) 

Of course, this invariance of the collision rules is a consequence of the instantaneous character 
of the hard collisions or, in other words, of the absence of an intrinsic time scale for hard 
particles. In the r time scale, the Boltzmann equation Q takes the form 

(1^ + ■ + ■ I:) = 

Here J is the same collision operator as in Eq. with but substituting the velocities v, Vi 
by w, Wi, and the scaled one-particle distribution function is 

fir,w,t) = iuJot)-'fir,v,t). (30) 

Therefore, the only modification of the Boltzmann equation is that a new term appears 
in the streaming part of the equation, as expected because of Eqs. fl27|l . This term has 
the same form as some thermostats introduced more or less artificially in order to allow the 
system to have a stationary state; however here it arises solely from a change in the time 
scale. 

Let us assume now that the system is initially in the HCS and that it remains in it along 
its time evolution. This implies to suppress the possibility that the system spontaneously 
develops spatial inhomogeneities due to the existence of a long wavelength hydrodynamic 



instability, the so-called clustering instability |15[. It is worth to stress that this effect 
is also present in the low density description provided by the Boltzmann equation 2^. 
In principle, the instability can be avoided by considering small enough systems, but this 
might lead to the presence of undesired finite size effects, especially for small values of 
the coefficient of restitution a. Nevertheless, at the level of description provided by the 
Boltzmann equation, this can be formally accomplished by restricting ourselves to consider 
the homogeneous form of Eq. (j^ . Then, a scaled temperature Tucsij) can be defined by 



d ~ f 1 ~ 

-nnkBTHcsiT) = J dw -mw^ fncsi'w, r), (31) 
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where fHcs{w,T) is the scaled one-particle distribution of the HCS. An evolution equation 
for this temperature is easily obtained by using Eq. (fTTj) . 



The solution of this equation is 

Tncsij) = 



I- - 2uJo] fncsir) = -(ff/^sir)- (32) 



c 



(33) 



It follows that any initial value of the temperature tends to a final steady value given by 

f.^(^J. (34) 

As discussed above, this result only holds as long as the system stays indefinitely in the 
HCS. By means of Eq. it can also be expressed in the equivalent form 

Co ^ (35) 

V0,st 

with 

vo,st = • (36) 

\ m J 

Therefore, ujo/vo^st is independent of Uq and proportional to an intrinsic property of the 
system, namely the dimensionless cooling rate Co- 

Moreover, once the scaled temperature has reached its steady value, the scaled distribu- 
tion function has the time independent form 

fst{w) = nnVo stXHCsic), c = = v. (37) 

Upon writing the last equality for c, we have taken into account the asymptotic form of 
Tucsif) given in Eq. (fT^ . The existence of a steady solution of Eq. (pUj) is enabled by the 
term proportional to cuq, so that the acceleration between collisions is able to balance the 
loss of energy in them. In order to avoid misunderstandings, it must be noted that the fact 
that the system is in the HCS does not impliy by itself that the temperature in the scaled 
variables takes its steady value. This only happens in the long time limit. However, what is 
relevant is that there is a mapping between the steady state in the scaled representation and 
the associated HCS, for arbitrary value of the parameter uoq. Of course, for any arbitrary 
state of the system, it is possible to relate every property in the original variables with the 
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corresponding (scaled) property in the scaled representation, but attention will be restricted 
in the following to the HCS. 

In order to confirm the practical usefulness of the steady representation, we have carried 
out DSMC simulations jlTI] of a low density granular gas whose dynamics is defined by 
Eqs. (|27|) and (|28|) . This A^-particle simulation algorithm is known to be consistent with 
the Boltzmann equation, in the sense that it provides numerical solutions of the equation. 
But the results coming from this kind of simulations go much further and cover all the 
properties of a dilute gas. This includes, in particular, fluctuations and non-equilibrium 
correlations, although the precise relationship between the simulation algorithm and the 
theory following from a description of the system in the context of kinetic theory, or non- 
equilibrium statistical mechanics, has not been established yet. In fact, the method tries 
to mimic, by means of a Markov process, the dynamics of a low density gas by uncoupling 
the streaming motion of the particles, given by Eqs. ()27|). and collisions during a small 
enough time interval, and also by neglecting velocity correlations in collisions. This is done 
independently of the number of particles being simulated (and, therefore, of their number 
density). In practice, this is implemented by dividing the coordinate space into cells of size 
smaller than the mean free path, and considering that all pairs of particles in the same cell 
can collide with a probability proportional to their relative velocity. 

One of the main technical advantages of DSMC, as compared with other particle sim- 
ulation methods, is that it allows to incorporate in the simulation algorithm the eventual 
symmetry properties of the particular physical situation of interest. The most trivial effect of 
this is the possibility of increasing the numerical statistics of the results, but it also permits 
to force the system to stay with a given symmetry, by eliminating from the dynamics the 
degrees of freedom associated with the "irrelevant" directions. For instance, for our present 
purposes we want the system to stay in the HCS, so that no spatial instabilities can be 
developed in any direction, and at the same time we want to avoid the introduction of finite 
size effects. This can be accomplished by simulating the A^-particle dynamics associated 
with the homogeneous Boltzmann equation. The way of implementing this in the simulation 
is by considering only one cell, i.e. any pair of particles in the system can collide and no 
attention is paid to their positions. Since the technical details of the method have been 



extensively discussed in the literature 







24| . they will not be reproduced here. 



The simulations we will present in the following correspond to a two-dimensional system 
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of hard disks, i.e. d = 2 . As already discussed, no boundary conditions are needed since we 
are simulating homogeneous situations and the positions of the particles play no role. For 
the same reason, there is no system size to be specified. Moreover, the numerical data we 
will report have been averaged over a number of trajectories, typically of the order of a few 
thousands. The number of particles in the system is = 10^, but we insist on that it has 
only a statistical meaning, since the dynamics of the system being used in the simulation 
corresponds always by construction to that of a low density gas. 

Accordingly with the scenario we have developed, the simulation of a low density gas 
whose underlying particle dynamics is defined by Eqs. (P7j) and (j2Hl), is expected to yield a 
steady state after a relatively short transient period. Averages of properties in the steady 
state are simply related with the (time dependent) properties of the HCS. In this way, the 
difficulties associated with the fast cooling of the fluid when described in terms of the actual 
variables, leading very soon to numerical inaccuracies, are overcome. One technical point 
requiring some attention is that the total momentum is unstable in the mapped represen- 
tation for all sizes of the system Round-off numerical errors lead to the presence of 
nonzero total momentum that grow very quickly due to the instability. This is easily avoid 
by calculating the total momentum at each simulation step and subtracting it evenly from 
the momentum of each particle. 

In the results to be reported in the following, the unit of time is given by 
£[2/cBT(0)/m]~^''^, the unit of length is i, the unit of mass is m, and has been set 
equal to unity. In Fig. [T] the evolution of the scaled temperature T is plotted as a function 
of the time r for several values of a, namely a = 0.5, 0.6, 0.7, and 0.9. The value of ujq 
used in each case is Uo = CT(0)^''^/2, with C{a) approximated by its estimate in the first 
Sonine approximation as given by Eqs. ()15|) and (j21I). If this latter expression were exact, 
the final steady temperature would be exactly the same as the initial one. The observed 
behavior is consistent with the theoretical predictions. In all the cases the temperature 
fiuctuates around a steady value after and initial transient time. In the original time scale, 
this corresponds to the regime in which the granular temperature has already reached its 
asymptotic form given in Eq. (fT^ . The steady value is very close to the initial temperature, 
as expected, although a small deviation is clearly identified for a = 0.5, indicating the ap- 
proximated nature of expression for (q. It is also seen in the figure that the amplitude of 
the long time temperature fiuctuations increases as the value of the coefficient of restitution 
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FIG. 1: Evolution of the scaled temperature T as a function of the scaled time r. Units are defined 
in the main text. The different line styles correspond to different values of a as indicated in the 
figure. 

decreases, i.e. as the system is more inelastic. This is due to the change of the shape of the 
velocity distribution function and also to the presence of velocity correlations in the HCS 
that are generated by the A^-particle dynamics of the system, even in the low density limit. 
This will be discussed in more detail elsewhere. 

/^From the values of Tgt, the cooling rate (o can be obtained as a function of a by means of 
Eq. ()35p. In Fig. |21these numerical results are compared with the theoretical prediction given 
by Eq. (j^ . As already indicated by the weak dependence on a of the steady temperature 
in Fig. there is a very good agreement. The discrepancy observed in the latter for a = 0.5 
can not be made out on the scale used in Fig. |21 

Since the steady velocity distribution has the same form as Xhcs: the steady represen- 
tation also provides very accurate data for it. As an example, in Fig. 01 the coefficient 02 
defined by Eq. (j7H) is plotted as a function of the coefficient of restitution and compared 
with the approximated expression given by Eq. ((221) • ^^is case, a small but systematic 
deviation is observed. The results presented so far in the above three figures are consistent, 
and physically equivalent, to those obtained previously by carrying out DSMC simulations 
in the actual variables of the system, i.e. by dealing directly with the time-dependent HCS 



22j |. The main advantage of the present representation is that it reduces the statistical 



uncertainties by introducing a steady state that maps exactly into the HCS. In the next 
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FIG. 2: The dimensionless cooling rate Co as a function of the coefficient of restitution a. The sym- 
bols are the values obtained from the simulations, while the solid line is the theoretical prediction 
in the first Sonine approximation. 

Section, it will be shown that the scaling is also useful to study time correlation functions. 

IV. AVERAGES AND TIME-CORRELATION FUNCTIONS IN THE STEADY 
REPRESENTATION 



Consider a phase function of the form 



N 



AiT)=J2aiRi,Vi), 



(38) 



i=l 



where F denotes a point in the phase space of the system and a(i?j, VJ) is a given one-particle 
dynamical variable. The average value of A at time t is 



{A;t) = J dT A{T,t)p{T,0) = J dV A{T)p{T,t). 



(39) 



Here p(F, t) is the A^-particle distribution function of the system at time t. This equation is 
equivalent to 

{A;t) = Jdr I dva{r,v)f{r,v,t), (40) 
that, using the scaling defined by Eqs. and (jSO)), can also be expressed as 



{A;t) 



dr / dw a 



r, (uot) ^w] f{r,w,T). 



(41) 
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FIG. 3: Coefficient 02 defined in Eq. (|2H) as a function of the coefficient of restitution a. The 
symbols are from the DSMC simulations, while the solid line is the theoretical prediction discussed 
in the main text. 



Suppose now a dilute system in the HCS. It has been established in Sec. Illll that in the 
long time limit, the scaled distribution function tends to the steady form fsti^) given by 
Eq. dnil), so that Eq. (gH) becomes 



{A; t)Hcs 
This is the low density limit of 

{A]t)HCS 

where f = {i?^, W^; z = 1, ■ ■ ■ , A^}, 



1 dr J dw a 


r, w 




Vo,st 



fst{w). 



(42) 



dV A{T)Pst{T), 



(43) 



N 



i=l 



vojt) 

Vo,st 



(44) 



and Pst(r) is the scaled A^-particle distribution of the steady state. 

In this way, averages in the HCS are exactly mapped onto averages in the the steady 
state of the scaled dynamics. It could be thought that since Eq. (j42j) has been obtained 
as the long time limit of Eq. ()41|1 for the HCS, its validity is somewhat restricted to very 
low actual temperatures Tncsit)- Nevertheless, this is not the case, because formally the 
initial temperature Thcs{^) can be as large as wanted and, consequently, Eq. can be 
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applied at any temperature Tucsif)- In fact, it is easily verified that the same equation 
can be derived from the particularization of Eq. ()40|) to the HCS by making the change of 
variable w = Vq^siC, without introducing any long time limit. Of course, the difference is 
that in this latter case THcsi^)-, present in Eq. (j42p through fo(t), can not be substituted 
by its asymptotic form. 

Next, let us consider time correlation functions in the HCS defined by 

CABit.t') = {A{t)B{ty,0)Hcs - {A;t)Hcs{B;t')Hcs, (45) 

with 

{A{t)B{t'y, o)hcs = jdv A(r, t)5(r, t')pHcs(r, o) (46) 

and it is assumed that t > t' > 0. The phase functions A{T) and -B(r) are the sum of 
one-particle dynamical variables, i.e. A is again given by Eq. (|38|l and 

N 

B{T) = J2b{R„Vi}. (47) 
Then, Eq. (j45p can be rewritten in the equivalent form 

CABit,t')= J dri J dvi J dr[ J dv[a{ri,Vi)b{r[,v[)hi/i{ri,Vi,t;r[,v[,t'), (48) 
where hi/i is the two-particle two-time correlation function of the HCS defined by Q| 

hi/iiri,vi,t;r[,v[,t') = fi/i{ri,Vi,t;r[,v[,t') - fHCs{vi,t)fHcs{v'i^t). (49) 
Here /i/i is the two-particle two-time reduced distribution function of the HCS, 

N N „ 
i=l j=l 

x5[r[-R,it)]6[v[-V,it)]. (50) 

In the low density limit, and by using the same kind of assumptions as needed to derive the 
Bofemaun equafo,, .t . possible to obtata a fonnal expressio,. for k,, involving o„ly the 
steady one-particle distribution function [16[. A sketch of the calculations is given in the 
Appendix |X1 When this expression is substituted into Eq. (PH|) it is obtained: 



CABit,t') = j dri J dwi fst{wi)a 



ri,- Wi,T - Ti 

V0,st 



ri,— wi 

Vo,st 



(51) 
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The time dependence of the dynamical variable a is given by 



a(ri,wi,r) = e^^^(-i)a(ri,^i^ 



(52) 



where Lb{wi) is the adjoint of the linearized Boltzmann operator around the steady state 
in the scaled dynamics, 



Lb(wi) = wi 



Ast = i^owi ■ 



d 
dri 



+ AstiWi) 



d 



dwi 



(53) 
(54) 



Kst{wi) = a'^ ^ j dw2fst{'W2) j dae{wi2-^)'Wi2-B[h„{wi,W2)-l]{l + Vi2). (55) 

The operator 7^12 interchanges the subindexes 1 and 2 to its right. Of course, if the dynamical 
variable a does not depend on the position r, Lb{wi) can be substituted by Kst{wi) in Eq. 

In the particular, but quite frequent case that a and h are homogeneous functions of the 
velocity of degree qi and g2, respectively, i.e. 



voit) 
ri,- — Wi 

Vo{t') 
ri, — Wi 



the time-correlation function becomes 



V0,st 
Vo,st 



a{ri,wi), 
b{ri,wi), 



CAB{t,t') = ° ^g^°g^ — / dri I dwi fst{wi)a{ri,Wi,T - Ti)b{ri,wi] 



(56) 



(57) 



In the context of the steady representation of the HCS of a low density granular gas as 
discussed in this paper, the relevant point of the above analysis is the following: suppose 
that a property of the gas can be expressed in terms of the time correlation function 



(a(r)6) 



St 



dr J dv fst{v)a{r,v,T)b{r,v) 



(58) 



with the time-dependence of a{r,v,T) given by Eq. (jK^ . A slight modification of the 
preceding discussion shows that this is the low density limit of 



where 



{A),t = J rfrp,,(r)A(r), (5),, = J rfrp,,(r)s(r). 



(59) 



(60) 
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{A{T)B;0),t = J dTpst{T)A{T,T)B{T), (61) 

with A{T) and -B(r) given by Eqs. and (jlTj), respectively. Besides, A(T,t) is generated 
from A{r) by the dynamics defined by Eqs. (j27jl and ()28|1 . and Pst(r) is the A^-particle 
distribution function of the steady state eventually reached by the system. Although the 
existence of this state is not rigorously proven, it is supported by the discussion carried in 
Sec, mil at the level of the Boltzmann equation and also by Molecular Dynamics simulations 
Q, I25I 1. The DSMC method provides an efficient tool to generate the A^-particle dynamical 
representation of a dilute gas, that is consistent with the Boltzmann equation. In summary, 
it allows the direct evaluation of Eq. (j59|) in the low density limit, where it is equivalent to 
Eq. (f5H|). 

An important application of the above ideas is the evaluation of the Navier-Stokes 
iransport coefficients of a dilute granular gas composed of hard spheres or disks. In refs. 
2 13; 27 1, these coefficients were derived from the inelastic Boltzmann equation by means 
of the Chapman- Enskog procedure, eigenfunctions expansions, and linear response theory, 
finding equivalent results. The expressions of all the transport coefficients are proportional 
to time integrals of correlation functions of the form 

Dab{s) = J dciXHCsici)a{ci,s)b{ci), (62) 

with the time dependence of a(c, s) determined by 

a(ci,s) = e^^^('^i)a(ci), (63) 

Ac(ci) = J dc2XHCs{c2) y"c/5- 0(ci2 ■ 5-)ci2 ■ 5-[6o.(ci,C2) - 1](1 + P12) + ^ci ■ (64) 
Finally, the time scale s is defined as 







'a'^. (65) 



Although the above representation is appropriate for formal manipulations, for com- 
putational purposes the dynamics associated with the operator Ac presents the technical 
complication that it involves the cooling rate (q, that therefore must be known a priori in- 
stead of being determined by the own simulation, as it is the case for the scaled dynamics 
being used in the present paper. For this reason, it is convenient to transform Eq. by 
writing it in terms of that dynamics. This is easily accomplished by introducing 

Co 2uoi ^ 

T = S, W = ——C = Vo,stC. (66) 
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Then, for an arbitrary function g{c) it is obtained 



w 



sK,g{c) = rAstg \^], (67) 
where A^t is the operator defined by Eq. ()54j) . It follows that Eq. (jH^ is the same as 

D^,{s) = Ua[^,^h[^),,, (68) 

N \V0,st J \V0,st/ 

where we use the notation introduced in Eq. (j58p . The above result relates the expression of 
the transport coefficients of a dilute granular gas, as obtained from the Boltzmann equation 
for the one-particle distribution function, with the low density limit of time-correlation func- 
tions in the A^-particle dynamics, computed in the steady state of the scaled representation 
introduced in Sec. IIIII 



V. SELF-DIFFUSION 



In ref . j28[ , the expression of the self-diffusion coefficient D of a dilute inelastic gas of hard 
particles has been derived from the Boltzmann-Lorentz equation by the Chapman-Enskog 



procedure. In Appendix ^ it is shown that the results reported in 
the form 

Vo{t)i 



28j can be expressed in 



D 



d 

voit) 



ds J dcxHcs{c)c{s) ■ c 

{w{t) ■ W)st. 



^ (69) 

Vo,stNd Jo 

In the last transformation we have used Eq. ()fi8|l . In fact, this is an special case, since the 
time dependence of c{s) is not given by the linearized Boltzmann collision operator Ac as 
in Eq. (jU^ . but by the Lorentz-Boltzmann one Asl,c(c), defined in Eq. ()B11|) . 



Consequently, in Eq. (jU^ it is 



cis 



W(T] 



AblA'^) = / dwi fstiwi) / daQ{g- B)g ■ B\h„{w, wi) - 1] 



Co d 
-w 



(70) 

(71) 

(72) 



2 ~~ dw 

Nevertheless, all the reasonings in Sec. |TV| and Appendix [X] can be easily adapted to the 
present case. The only, but relevant, difference, is that now the expression {w{t) ■ w)st is 
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FIG. 4: Time evolution of the normalized velocity autocorrelation function Cto(t)/Cot(0) as a 
function of the scaled time r. The results for three different values of the coefficient of restitution 
a are shown, as indicated in the figure. 

the low density limit of 

N 

a,(r) = J2l dV pstiT)W,it) ■ W, (73) 
1=1 

and not of 

N N „ 

d^Psti^)W^{t)■W,. (74) 

i=l j=l 

Of course, this has to be taken into account when computing the self-diffusion coefficient by 
means of DSMC simulations. 

Figure m shows the velocity autocorrelation function C^,^,(r) for three different values of the 
coefficient of restitution. It is seen that it decays quite fast to zero, implying the convergence 
of its time integral and, therefore, the existence of the constant diffusion coefficient. The 
value of the velocity autocorrelation at each time r is based on an average of all the times 
included in the simulation, once the steady state has been reached. In Fig. the same 
quantity is plotted on a logarithmic scale where an exponential decay for all times is clearly 
identified. 

/^From the velocity autocorrelation function the self-diffusion coefficient is obtained by 
means of the Green-Kubo relation ()69|) . Figure IHl shows the ratio of the obtained values 
D{Thcs,C() to the elastic limit (a 1) of the value predicted by the Chapman- Enskog 
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X 



FIG. 5: The same as in Fig. |l]but on a logarithmic scale. 

method in the first Sonine approximation Dq{Thcs)i namely 

D (a) = -jrjTf, ^, (75) 

T{d/2)d ( kBTHcs \"' 
MThcs) = [-^) ■ (76) 

Also plotted is the theoretical prediction for D*{q) again in the first Sonine approximation 



D*(a) 



;i + «f-^(4 + a-3a^) 
lo 



-1 



(77) 



The agreement between theory and simulation is quite good in all the range of values of 
a considered, although a systematic deviation, which increases as the value of a decreases, 
is observed. Self-diffusion in the HCS of a system of inelastic hard spheres {d = 3) has been 
previously studied by means of DSMC in the actual variables, i.e. under continuous cooling, 
in ref . j28| , where the diffusion constant was calculated from the mean-square displacement 
by means of the (inelastic) Einstein relation and also from the linear response of the system 
to a density perturbation. The comparison between the numerical results obtained there 
and the theoretical prediction given by Eq. (j77|) is very similar to the one presented in Fig. 
EJ Very recently ^2^, it has been shown that the agreement improves over the whole range of 
values of a if D*{a) is computed in the second Sonine approximation. Although the analysis 
is restricted to the case of inelastic hard spheres, it seems sensible to expect the same kind 
of behavior for a system of inelastic hard disks. 
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FIG. 6: Dimensionless reduced diffusion coefficient D* as a function of the coefficient of restitution 
a. The symbols are from the simulations, while the solid line is the theoretical prediction from the 
Chapman-Enskog procedure in the first Sonine approximation. 

As already mentioned several times, the main feature of the numerical method used here 
is that it takes advantage of the steady representation of the HCS, removing any limitation 
on the time on which trajectories of the system may be followed. In addition, let us point 
out that the analysis of the self-diffusion coefficient based on the velocity autocorrelation 
function as developed here, is expected to give more accurate results, from a statistical point 
of view, than those based on the mean square displacement. This is because, while each 
value of the former in a given trajectory of the system is obtained from an average over all 
the simulation interval, each value of the latter is obtained from a single evaluation. 

VI. DISCUSSION AND CONCLUSION 

In this paper, the actual dynamics of a low density granular gas composed of smooth 
hard inelastic particles has been transformed to another one in which the HCS becomes 
time-independent. In this way, the need for more or less uncontrolled mechanisms such as 
internal thermostats in order to get a stationary state is eliminated. The transformation is 
closely related to the fact that the temperature of the HCS becomes independent from its 
initial value in the asymptotically long time limit, a property that has not received enough 
attention up to now. The exact correspondence between both formulations for averages 
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and time correlation functions has been explicitly established. This requires to consider 
the kinetic equation for the one-particle distribution function, i.e. the inelastic Boltzmann 
equation, and also the equation for the two-particle and two-time correlation function in 
the low density limit. The latter is obtained by extending in a natural way the standard 
methods of kinetic theory. In particular, it has been shown that the steady temperature 
directly determines the value of the cooling rate of the system. Then, the DSMC method has 
been used to measure it and the numerical results have been compared with the theoretical 
predictions obtained by solving the Boltzmann equation in the first Sonine approximation. 

The introduction of the steady representation of the HCS enables the evaluation of the 
Green-Kubo expressions for the transport coefficients of a dilute granular gas by means of 
the DSMC method, just as for normal fluids whose particles collide elastically. To put this 
in a proper context, we have emphasized that the DSMC method is not just a numerical 
tool to solve the Boltzmann equation, but it formulates an effective N-body dynamics that 
is expected to be equivalent to the Newtonian one in the low density limit. Moreover, it 
has the advantage of allowing to incorporate in the own effective dynamics of the parti- 
cles the symmetry properties of the particular state being simulated. So, it is possible to 
force the system to stay homogeneous, avoiding the spontaneous development of spacial 
inhomogeneities associated with the clustering instability. Here this has been illustrated 
for the simplest case of self-diffusion, whose expression involves the velocity autocorrela- 
tion function. The results have been compared with the theoretical prediction obtained 
by the Chapman- Enskog method and also with previous numerical simulations carried out 
in the actual dynamics in which the HCS is time-dependent. The study of the remaining 
Navier-Stokes transport coefficients will be reported elsewhere. 

In summary, we believe that the transformed dynamics discussed in this work is of both 
formal and practical interest for the study of granular fluids in the low density limit. Al- 
though we have restricted here ourselves to its application to the HCS, the mapping with 
the actual dynamics can be easily extended to any arbitrary state. 
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APPENDIX A: LOW DENSITY LIMIT OF TIME-CORRELATION FUNC- 
TIONS IN THE HCS 



In the low density limit, neglecting three-particle correlations, the function hi/i obeys 
the equation, 

hi/i{ri,vi,t]r[,v[,t') = (Al) 



where K is the linearized Bolztmann operator 

K{vi,t) = a^-^ j dv2 I daeivi2 ■ ^)\vi2 ■ a\ia-%„ ~ +Vu)fHcs{v2,t). (A2) 

The permutation operator Vu interchanges the indexes 1 and 2 of the velocities appearing 
to its right. The above equation holds for t > t' > and it has to be solved with the initial 
condition 

hi/i{ri,vi,t';r[,v[,t') = g2,Hcs{ruVi,r[,v[,t') + 6{ri - r[) 6{vi - v[)fHcs{vut'), (A3) 

g2,Hcs being the two-particle one-time correlation function of the HCS. Equation (jAljl can 
be derived by the hierarchy method starting from the pseudo-Liouville equation for a system 
of inelastic hard spheres or disks. A detailed analysis for the elastic case is presented in ref. 
[igI i. Since the same method can be applied mutatis mutandis to the inelastic case, it will 
be not reproduced here. Now, time and velocities are scaled by 



and 



'0 Vo,st 



Vo,st 
W = — V, 



vo{t) 



(A4) 



(A5) 



respectively. We are using the same symbols as in Sec. IIIII since we have seen that in the 
long time limit the scaling defined by Eq. ()A4|1 is equivalent to that defined by Eq. 
The scaled correlation function is 

, 1 d 



hi/i{ri,Vi,t;r[,v[,t') 



(A6) 



and, in terms of the new variables, Eq. (jAlj) becomes 

hi/i{ri, wi, r; r[, w[, t') = 0, 



9 d . . . 

— + wi- Ast{wi^ 

or or I 



(A7) 
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valid for t > t'. The operator Agt is defined by 



d 



Astiwi) = Kst{wi) - ujQ- it>i, 

OWi 

where Kgt is the steady hnear Boltzmann coUision operator 



(A8) 



Kst{wi) = a''-^ j dw2 I daQ{wi2-^)wi2-a[a~%\wi,W2)-l\ (l + Vi2)fstiw2). (A9) 



Integration of Eq. ()A7|) gives 

/ii/i(ri, w,, r; r[, w[, r') = e(^-^')^^('"^)/ii/i(ri, w^, r'; r[, w[, r'), 



Lb{wi) = Astiwi) - wi 
The initial condition on the right hand side is 



d 

dvi 



(AlO) 
(All) 



Wi, r'; r[, w[, r) = g2,HCs{ri, Wi,r[, w[, r') +5(^1 - r[) 6{wi-w[)fstiwi), (A12) 



with 



92,Hcs{ru'Wi,r[,w[,T') 



vojt') 



-I 2d 



92iri,vi,r[,v[,t'). 



(A13) 



Now the assumption is made that the contribution to Eq. (jAlOj) coming from g2 can be 
neglected in the low density limit we are considering. Although there is no explicit proof 
for this, it is consistent with the hypothesis made to derive the Boltzmann equation. On 
the other hand, it must be realized that the HCS is not an equilibrium state and, therefore, 
it can present relevant position and velocity correlations, but we expect them to become 
negligible for asymptotically small densities. Then, substitution of Eq. ()A12|) into Eq. 
yields 



CAsit^t') = j dri j dwi J d.r[ J dw[a 



voit) 
ri,— — wi 

Vo.st 





Vo,st 



dvi / dwi / dr[ / dw[ b 



V0,st 



5{ri - r[) 5{wi - w[)fst{wi] 



g(r-r')Ls{^l)^ 



ri,- — wi 



V0,st 

dri I dwifst{wi)b 



n, -3 — wi 

V0,st 



^(T~T')LB{w-i_)^ 



Vo{t) 
ri,- Wi 



(A14) 



where is the adjoint of and it is given by Eq. The above equation is the same 

as Eq. dni- 



25 



APPENDIX B: GREEN-KUBO EXPRESSION FOR SELF-DIFFUSION FROM 
THE CHAPMAN-ENSKOG RESULTS 



In ref. j28[, the self-diffusion coefficient of a dilute granular fluid was obtained from 
the Boltzmann-Lorentz by the Chapmann-Enskog method. Equations (24) and (26) in the 
above-mentioned reference are 



D = ~- f dvv- B(v), 
a J 



where B{v) is the solution of the integral equation 

d 



Kbl + ChcsThcs 



1 



B{v) = — vfncsiv) 
oThcs J riH 



We introduce dimensionless time and velocity scales by 



s = I dt' c 

I 



V 



In terms of them, Eq. ()B2|) becomes 



Kbl,c - ■ cj B{c) = cxHcsic), 



with 



and 



B(c) 



Kbl Ac) 



vt{t) 



vo{t) 



B(v) 



Kbl{v). 



The formal solution of Eq. ()B5j) can be written as 

roo 

Jo 

^blA^) = KblA^) - ■ 

2 oc 



Substitution of Eq. ()B8|1 into Eq. ()B1|1 gives 



d 



d Jo 



ds / dc 



■ cxhcs[c) 



iBl) 



(B2) 



Kbl being the inelastic Boltzmann-Lorentz collision operator 

Kbl{v) = cr"-' I dvi f dae{g ■ a)g ■ a \a-%-\v,Vi) - l] fHcs{vi,t). (B3) 



(B4) 

(B5) 

(B6) 
(B7) 

(B8) 
(B9) 



(BIO) 
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where h-BL,c{c) is the adjoint of h-BL,c{c), 

Abl,c(ci) = j dc2 Xhcs{c2) j da- e(ci2 • o-)ci2 • a-[6o-(ci, C2) - 1] + ^Ci • . (Bll) 
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